function Mlbar_sim = Sim_Firmsize(pm,p,w,ze,zr)

    %If the guess converges, then we proceed with the data moments
    %Average x
    sim.Lbar_i = zeros(pm.I,1);
    sim.Lbar_f = zeros(pm.I,1);
    
    for r = 1:pm.I
        sim.Lbar_i(r) = integral(@(e) integral(@(x)  l_i(pm,x,e,p(r),w,r)...
                .*lognpdf(x,0,pm.sigma).*lognpdf(e,0,pm.se),ze(r),zr(r),'ArrayValued',true),-Inf,Inf) ...
                ./(logncdf(zr(r),0,pm.sigma)-logncdf(ze(r),0,pm.sigma));

            sim.Lbar_f(r) = integral(@(e) integral(@(x)  l_f(pm,x,e,p(r),w,pm.t(r),r)...
                .*lognpdf(x,0,pm.sigma).*lognpdf(e,0,pm.se),zr(r),Inf,'ArrayValued',true),-Inf,Inf) ...
                ./(1 - logncdf(zr(r),0,pm.sigma));
    end
    %BC and BP for the industries
    %    sim.lbar_i = 1./(pm.sigma-pm.rho).*(log(pm.rho./pm.sigma.*p./w) + zbar_i);
    %    sim.lbar_f = 1./(1-pm.rho).*(log(pm.rho.*(1-pm.t).*p./wf) + zbar_f);
    Mlbar_sim = [sim.Lbar_i;sim.Lbar_f];    
